import numpy as np
import matplotlib.pyplot as plt

# 有一组血糖指数X与血脂指数Y之间的测试数据如下:
X = [4, 3, 3, 4, 2, 2, 0, 1, 2, 5, 1, 2, 5, 1, 3]
Y = [8, 6, 6, 7, 4, 4, 2, 4, 5, 9, 3, 4, 8, 3, 6]
x = np.array(X, dtype=np.float64)
y = np.array(Y, dtype=np.float64)
m = len(x)

# 放缩：标准化
xy = np.c_[x, y]
mu = xy.mean()
sigma = xy.std()
xy -= mu
xy /= sigma


def model(x, theta):
    """
    Get hypothesis vector

    :param x: matrix of x
    :param theta: theta vector
    :return: hypothesis vector
    """
    return x.dot(theta)

#
# 请通过Python实现一元线性回归模型，并用此模型预测y，具体要求如下：
#
# 实现线性回归的代价函数
def cost_function(h, y):
    """
    Get cost function value

    :param h: hypothesis value vector
    :param y: target value vector
    :return: cost function value scalar
    """
    m = len(h)
    e = h - y  # mx1
    sq = e.T.dot(e)  # 1x1
    sq = sq[0][0]
    cost = 1.0 / (2.0 * m) * sq  # scalar
    return cost


# 实现梯度下降函数
# 要求输出迭代过程中的代价函数值
def gradient_descent(x, y, theta=None, alpha=0.001, num_iter=1500):
    """
    Gradient descent algorithm

    :param x: matrix of x
    :param y: vector of target values
    :param theta: theta vector
    :param alpha:
    :param num_iter:
    :return:
    """
    m = x.shape[0]
    n = x.shape[1]
    if theta is None:
        theta = np.zeros([n, 1])

    j_hist = []
    for i in range(num_iter):
        h = model(x, theta)
        e = h - y
        j = cost_function(h, y)
        j_hist.append(j)

        delta = alpha * (1.0 / m * x.T.dot(e))
        if np.allclose(np.zeros_like(delta), delta, atol=1e-3):
            return theta, j_hist, h.ravel()

        theta -= delta
    return theta, j_hist, h.ravel()


# 通过梯度下降计算回归模型，用所得模型对测试集的数据进行预测
num_iter_ori = num_iter = 1500
theta, j_hist, h = gradient_descent(x.reshape(m, 1), y.reshape(m, 1), num_iter=num_iter)
if num_iter_ori > num_iter:
    print(f'Converged at {num_iter}th iteration.')
else:
    print(f'Not converged after {num_iter} iterations.')
print(f'theta = {theta}')
plt.figure(figsize=[12, 5])
plt.subplot(1, 2, 1)
plt.plot(j_hist)
plt.grid()

# 以横轴为真实值，纵轴为预测值，画出散点图进行对比
plt.subplot(1, 2, 2)
plt.scatter(y, y, c='b', label='Real value')
plt.scatter(y, h, c='r', label='Hypothesis')
plt.legend()
plt.grid()

plt.show()
